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We develop an effective mass theory for substitutional donors in silicon in an inhomogeneous 
environment. Valley-orbit coupling is included perturbatively. We specifically consider the Stark 
effect in Si:P. In this case, the theory becomes more accurate at high fields, while it is designed to 
give correct experimental binding energies at zero field. Unexpectedly, the ground state energy for 
the donor electron is found to increase with electric field as a consequence of spectrum narrowing 
of the Is manifold. Our results are of particular importance for the Kane quantum computer. 



Phosphorus doped silicon is one of the most well- 
studied semiconducting systems Q, owing to its impor- 
tance in the electronics industry. A more recent and 
exotic application is quantum computing, in which the 
Kane proposal posits nuclear spin qubits on Si:P, with 
interactions mediated by donor-bound electrons [2(. By 
tuning the potentials on aligned electrodes, the electrons 
are brought to the point of ionization, in order to control 
nuclear hyperfine interactions and the overlap between 
neighboring qubits. Spin dependent ionization also pro- 
vides a means for electrical detection of the qubit state. 
The precise control of donor-bound electrons in such a 
complex environment remains an experimental challenge 
Q , requiring detailed theoretical input 0, 0, @ ■ 

The theory of isolated donors in silicon remains one of 
the crowning achievements of solid-state physics. One of 
the most effective treatments for shallow donors is the 
effective mass approximation (EMA) Q, in which the 
donor potential is assumed to vary slowly compared to 
the crystal potential. As a result, the long and short- 
wavelength physics decouple. An excess electron on the 
donor can be described by an envelope equation - a 
Schrodinger equation with an effective mass and a dielec- 
tric constant. The theory highlights the most essential 
feature of silicon's bulk band structure: the six-fold de- 
generacy of the conduction valleys. Near the impurity 
core, the assumption of slowly varying potential breaks 
down, causing valley-orbit interactions to couple enve- 
lope functions in different valleys. A careful treatment 
of the potential very near the impurity (the 'central cell' 
region) enables estimates for the energy splitting of the 
six valley states, which are in good agreement with ex- 
periments ||- 

In a more general, inhomogeneous environment, sili- 
con donors can be studied by tight-binding techniques 
0. However, existing EMA theories introduce severe 
approximations that limit their scope. Many important 
questions therefore remain open. For example, there is no 
theory to determine how the weight of a donor wavefunc- 
tion will redistribute among the six conduction valleys in 
the presence of a generic (i.e., low symmetry) potential. 
Since the envelope functions in different valleys have dif- 
ferent energies (due to anisotropy of the effective mass), 
this is an important question. 



In the context of quantum computing, several recent 
papers obtain tractible results for donor ionization, by 
ignoring valley-orbit interactions 0, 0, El- This single- 
valley approach provides a useful picture of the distorted 
envelopes in the individual valleys. However, it does not 
capture the spectrum narrowing of the Is manifold. Smit 
et al. have recently studied the Stark effect, beyond the 
single- valley EMA, byapplying symmetry arguments and 
perturbation theory |lfj |. However, their results are only 
applicable at low fields, and they also do not capture 
spectrum narrowing. In the present Letter, we show how 
to overcome such difficulties. We develop a multi-valley 
EMA for shallow donors in silicon in a general, inhomo- 
geneous environment. 

Before outlining the theory, we discuss two compet- 
ing effects that determine the energy shift of electrons 
in the Is manifold. The most well known effect is the 
quadratic Stark shift, which causes the average energy of 
the six valleys to decrease with field. The second effect is 
spectrum narrowing within the Is manifold. Valley-orbit 
effects induce a manifold splitting of about 13 meV for 
the six Is states in Si:P. As the field increases, the donor 
electron is gradually pulled off the impurity and out of 
the central cell, causing the splitting to narrow. Con- 
sequently, the ground state shifts upwards toward the 
average energy of the manifold. An accurate, quantita- 
tive analysis is needed to establish the dominant effect. 
The energies computed in our EMA analysis are shown 
in Fig. ^ We find that spectrum narrowing is the dom- 
inant effect for the ground state, causing the energy to 
increase with field. A similar result was obtained re- 
cently in tight binding simulations of a gated qubit 0]. 
However, the differences in that geometry make compari- 
son difficult. Here, the surprising behavior arises directly 
from the multi-valley band structure. 

We now describe the effective mass theory. The wave- 
function of the excess electron of the donor impurity can 
be written as 
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*(r)=2a i i (r)fl(r) 1 (1) 
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where at are the valley composition parameters, reflect- 
ing the portion of the wavefunction in each of the six 
valleys. Normalization gives 
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functions are given by 4>i{r) — Ui(r)e lki ' T , where k; = k$i 
specifies the minimum of the ith conduction valley, and 
Ui has the periodicity of the crystal lattice. The Fi are 
envelope functions for the six valleys. 

To develop an envelope equation for Fi, we consider the 
potential energy U(r), including both the impurity Vl(r) 
and its surroundings, but excluding the crystal poten- 
tial Q. In this Letter, we specifically consider a uniform 
external field E = Ez. However, the results are eas- 
ily generalized. The impurity ion is not a perfect point 
charge. Deviations from point charge behavior U cc (r) are 
called 'central cell corrections,' because they are strongly 
localized within a central cell radius of 1-2 A ■ Thus, 
U(r) = Vi(r) - eEz = [-e 2 /4ner + U cc (r)] - eEz. Cen- 
tral cell effects can also include the breakdown of the 
concept of the dielectric constant and local distortions of 
the Si lattice near the P ion. Although there has been 
progress in the analytical description of the central cell 
|l2| ) a full elucidation is challenging. Below, we overcome 
this difficulty by introducing a semi-analytical treatment 
of central cell effects, using the known energy spectrum 
of the Si:P donor electron at zero field. 

The envelope equation of Fritzsche and Twose can be 
extended to include an external field. The equation goes 
beyond the sing le- valley EMA by including valley-orbit 
interactions [I1E1E1E3: 
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FIG. 1: Solid curves: the variational binding energy as a func- 
tion of electric field for the Is manifold in Si:P. (The second 
solid curve from bottom is doubly degenerate.) Dashed curve: 
average energy of the manifold. 



perturbation, we solve the uncoupled, single-valley equa- 
tions to obtain the eigenfunctions Fj(r). These are 
then used to compute the valley-orbit terms. The result- 
ing first-order Hamiltonian is given in the valley basis 
(+x,-x,+y,-y,+z,-z): 



j'=i 



Aner 



eEz-e}Fj(r) 



+ J2a j e ik f r U cc (r)F j (r) = 0. (2) 
i=i 

We refer to the first line of J2J as EMA terms, and the 
second line as central cell corrections. Because the central 
cell terms are localized so strongly near the impurity, 
U cc (r) is essentially a contact potential 0|- This allows 
us to replace Fj(r) by Fqj, the envelope amplitude at the 
impurity site. 

In the so-called multi- valley EMA of Eq. @ , e is the 
energy eigenvalue. The kinetic energy operator Tj(fe) is 
the quadratic expansion of the conduction band disper- 
sion relation E c (k) around the valley minimum fcj. This 
prescription for Tj introduces spurious inter-valley con- 
tributions, which can be avoided if necessary 0, Il7| . 
However, in our theory, only the intra-valley energy ex- 
pression will be computed explicitly For this case, the 
quadratic expansion is appropriate. Note that we have 
limited our consideration to the six low-lyingvalleys in 
the Is manifold, as appropriate for the EMA Higher- 
bands give contributions much smaller than the terms 
considered here. 

To solve the envelope equation we take a perturbative 
approach by assuming that valley-orbit interactions are 
weak. We shall justify this assumption later, by com- 
paring the sizes of various terms. At zeroth order in the 
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where E% = E t + A oi , and we have used axial (E x = E y ) 
and time reversal (E +z = E- z ) symmetry. The various 
terms in H are defined as follows: 
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Equation JSJl also generates off-diagonal terms in the 
Hamiltonian, analogous to Because of the fast oscil- 
lations in the off-diagonal integrals and the smoothness 
of the effective mass eigenfunctions at zeroth order, such 
terms can be ignored in comparison with Ej The 
only non- vanishing off-diagonal terms are Ay and A2 ji, 
due to the strong localization of U cc on atomic length 
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scales. Approximating U cc as a contact potential (l5| 
gives 

A 0j = v F 2 p Ay = Vl F^, A 2ji = v 2 F 0j F 0h (8) 

where Foj is the magnitude of the envelope function at 
the impurity site. Note that we have chosen to include 
the electric field in the Ej terms, rather than the central 
cell terms. This assignment is justified by comparing 
—eEz with the point charge potential — e 2 /Aner. Within 
the central cell, the former becomes dominant only for 
fields much larger than the ionization field. The contact 
potentials i>o,i,2 are therefore independent of field. They 
can be determined empirically at zero field, where the 
energy spectrum of the donor is well known. This is a 
central observation of our theory. Once i>o,i,2 are speci- 
fied, the field dependence of the valley-orbit terms arises 
only from the wavefunction normalization Fq? . Since 

Fq^ decreases with increasing field, valley splitting also 
decreases. This is the origin of spectrum narrowing. 

We now perform the Stark shift calculation for Si:P 
using the pertubation theory. The unperturbed (single- 
valley) envelope equations are given by 
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The equation for Fy(r) is identical to F^ (r), with 
x «-> y. The longitudinal and transverse effective masses 
for silicon are rn[ — 0.916m and m* t = 0.191tooi re ~ 
spectively. The low temperature dielectric constant is 
e = 11.4e . 

We use a variational method to solve Eqs. ® and 1)11)0. 
obtaining good agreement with numerical results. The 
variational forms are given by 

Fj°)(r) = F 0x (l + q x z) exp (- V*7^ + iv 2 + ^ 2 )/^) , 

F^(r) = F 0z (l + q z z) exp (-^ z 2 / a 2 + {x 2 + y 2 ) /b 2 ) , 

These functions are similar to ones used successfully in 
the zero- field case At this order, the envelope equa- 
tions are uncoupled, and the valley energies E x °l are min- 
imized independently. 

The contact potentials are determined at zero field. In 
this limit, the problem is isotropic, giving E x ^ — E z ^ = 
—31.28 meV, a x = a z = 1.360 nm, b x = b z = 2.365 nm, 
q x = q z = 0, and F 0x = F 0z = 6.469 x 10 12 m" 3 / 2 . 
The central cell terms do not require valley indices: 
Aq,Ai,A2. The zero- field binding energies are well 



known from experiments Q: e = —45.59(1), —33.89(3), 
and —32.58(2) meV (parentheses indicate the zero-field 
degeneracies). By requiring the Hamiltonian @ to pro- 
vide these eigenvalues, the matrix elements are uniquely 
specified: E = -35.40 meV, A = -4.13 meV, Ai = 
— 1.51 meV, A2 = —2.17 meV, leading to the identifica- 
tion v = -1.58 x 10~ 47 Jm 3 , vi = -5.78 x 10~ 48 Jm 3 , 
and v 2 = -8.31 x 10~ 48 Jm 3 , from Eq. ®. Note that 
these values depend on our choice of variational func- 
tions. We can now justify the perturbation approach by 
comparing the perturbation terms A04.2 to the EMA en- 
ergy £?(°). The theory improves in accuracy with increas- 
ing field, as the central cell terms grow smaller. However, 
the method is also designed to obtain exact experimental 
results for the donor binding energies at zero-field. 

At non-zero fields, we use the same perturbation pre- 
scription. The variational form is used to minimize the 
energies of the uncoupled envelope functions, g iving E x \ 
and F x ® z . The off-diagonal elements of H are obtained 
from JHJ, using the field-independent values of vo.1,2 just 
obtained. 

We diagonalize H to obtain the first-order eigenstates. 
The resulting eigenvectors are the valley composition pa- 
rameters ai. Using the zero-field notation of Ref. [lOj, 
we obtain the lowest (g) and highest (r) eigenvalues: 



1 r ~ 



E x + E z + Ai x + Aiz + 2A 2£! ;^ 



(11) 



±J32A 2 XZ + {E x ~E Z + A lx - A lz + 2A 2xy ) 2 . 



Corresponding eigenvectors are expressed in the valley 
basis: 



(l, 1, 1, 1, a' gi r, o / 9 ,r)/y / 4 + 2a' 
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9,r ; 



(12) 



where 



9 > r ~ 4A, 



(E x -E z + A lx - A lz + 2A 2xy ) 



±^32A 2 XZ + {E x -E z + A lx - A lz + 2A 2xy ) 2 

(13) 

At zero field, we recover the expected eigenstates a g = 
(1, 1, 1, 1, 1, 1)/V6 and OL r = (1, 1, 1, 1, -2, -2)/yT2. 

The remaining eigenvectors are the same as their zero- 
field expressions: 



OL x 


= (1,-1,0,0,0,0)/^, 
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= (0,0,1,-1,0,0)/V2, 


(15) 


OL z 
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with eigenvalues, from lowest to highest: 

£x,y — E x — Ai x , (18) 
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FIG. 2: Amplitude of the ground state envelope function vs. 
electric field. Inset: ground state valley composition parame- 
ters. The solid lines correspond to the x valley. The dashed 
lines correspond to the z valley. All curves are normalized to 
1 at zero field. 



potential. The theory is applied to the Stark shift in Si:P. 
In contrast with previous theories, we predict that the 
ground state energy of the donor will increase with elec- 
tric field, due to spectrum narrowing of the Is manifold. 
Comparison with previous results [10| gives corrections 
as large as 10 meV. Such theories also do not address 
the valley redistribution of the wavefunction, which can 
amount to 15%. The new effects are most prevalent in 
the ionization regime, where quantum computers are ex- 
pected to operate. These remarkable results show that 
even qualitative descriptions of shallow donors in Si must 
take into account valley-orbit interactions. 
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Foundation through the ITR program (DMR-0325634) 
and the QuBIC program (EIA-0130400). 



e z =E z -A lz , (19) 

e s = E x + A lx - 2A 2xy . (20) 

The resulting energies, obtained from this 
perturbation-variational theory, are plotted in Fig. ^ 
For increasing fields, the wavefunction moves off the 
impurity site, as shown in Fig. [21 In both figures, the 
results are displayed for fields up to the critical field 
(~ 3.7 MV/m), beyond which the system is completely 
unstable to ionization, and solutions cannot be obtained. 
Because of the observed upturn of the ground state en- 
ergy with electric field, ionization occurs at considerably 
lower fields than expected from a single-valley EMA. 
The latter predicts a downturn of e g with field. Clearly 
multi-valley effects play an crucial role in ionization 
calculations. 

An interesting question is whether the electric field can 
redistribute the weight of the electron between the six 
valleys. As shown in the inset of Fig. |3 the answer is 
"yes." The effect is small, except near the critical field. 
However, for many donor-bound qubit schemes, the ion- 
ized or nearly-ionized donor states are utilized for gate 
operations 0, Hil H^ |. so the regime is of practical in- 
terest. Valley redistribution is also of interest when an 
unintended impurity interacts with a quantum dot qubit 
in a quantum well [2jj. Because of strain effects, the 
qubit wavefunction comprises only the two z valleys [2l| , 
so an exchange coupling between the qubit and impurity 
electrons occurs primarily in the z valley channel. It is 
therefore necessary to know what portion of the donor 
wavefunction resides in the z valleys. 

In conclusion, we have developed a multi-valley EMA 
for donor-bound electrons in silicon in an inhomogeneous 
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